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Abstract 


An  unstructured  grid  Large  Eddy  Simulation  (LES)  methodology  has  been  developed  for  compress¬ 
ible  high  speed  flows.  The  filtered  compressible  Navier-Stokes  equations  are  solved  on  an  unstruc¬ 
tured  grid  of  tetrahedra.  The  inviscid  fluxes  are  obtained  from  an  exact  locally  one-dimensional 
Riemann  solver  using  Godunov’s  method.  The  viscous  fluxes  are  obtained  using  a  discrete  analog  of 
Gauss’  Theorem.  The  reconstruction  is  performed  using  a  Least  Squares  technique.  The  temporal 
integration  is  a  Runge-Kutta  method.  The  algorithm  is  overall  second  order  accurate  in  space  and 
time.  Four  flowfields  have  been  computed:  supersonic  flat  plate  boundary  layer,  supersonic  com¬ 
pression  corner,  supersonic  expansion-compression  corner  and  subsonic  square  jet.  The  computed 
results  show  close  agreement  with  experiment  and  Direct  Numerical  Simulation,  and  validate  the 
unstructured  grid  LES  methodology. 
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Introduction 


The  effective  design  of  high  speed  aircraft  and  missiles  depends  critically  upon  accurate  prediction 
of  aerodynamic  and  aerothermodynamic  performance  which  are  strongly  affected  by  flow  turbulence 
under  most  flight  conditions.  Prom  an  engineering  standpoint,  the  aircraft  or  missile  aerodynamicist 
needs  the  capability  for  accurate  prediction  of  the  mean  and  rms  fluctuating  surface  pressure  (pw 
and  p'w)  and  surface  heat  transfer  (qw  and  q'w),  mean  surface  skin  friction  (tw),  and  locations  of 
primary  and  secondary  separation. 


Table  1:  RANS  Capability  for  3-D  Shock  Wave  Boundary  Layer  Interaction 


Quantity 

Satisfactory 

Unsatisfactory 

No  capability  shown 

Pw 

V 

p'w 

<lw 

V 

Qw 

Primary  Separation 

V 

Secondary  Separation 

V 

The  current  methodology  for  prediction  of  compressible  turbulent  flows  is  based  on  the  Reynolds- 
averaged  Navier-Stokes  (RANS)  equations  (Knight  1993).  This  approach  has  yielded  a  hierarchy 
of  turbulence  models  extending  from  zero-equation  to  full  Reynolds  Stress  Equation  models.  While 
these  models  have  generally  been  capable  of  predicting  the  engineering  quantities  of  interest  in 
weakly  perturbed  boundary  layers,  they  have  been  unable  to  accurately  predict  the  complex  3-D 
flows  which  are  encountered  in  highly  maneuvering,  high  angle-of-attack  flight.  Two  recent  extensive 
reviews  have  documented  the  capabilities  and  deficiencies  of  a  wide  range  of  RANS  models  for 
prediction  of  complex  3-D  flows  with  shock  wave-turbulent  boundary  layer  interactions  (Knight 
1997,  Knight  and  Degrez  1998).  The  results,  summarized  in  Table  1,  indicate  that  a  significant 
number  of  critical  engineering  quantities  are  not  capable  of  prediction  by  current  RANS  models. 
Therefore,  more  advanced  turbulence  models  are  needed  which  have  the  ability  to  simulate  the 
complex  physics  of  turbulence  with  greater  generality. 

Large  Eddy  Simulation  (LES)  is  an  alternative  to  RANS  which  may  be  capable  of  predicting  more 
(or  all)  of  the  aerodynamic  and  aerothermodynamic  quantities  of  engineering  interest  described 
above.  In  LES,  the  governing  equations  are  spatially  filtered  on  the  scale  of  the  numerical  grid.  The 
large,  energy-containing  eddies  are  directly  computed.  These  eddies  are  strongly  influenced  by  the 
physical  geometry  and  configuration  of  the  flow.  Thus,  the  direct  computation  of  the  large  eddies  by 
LES,  as  opposed  to  the  modeling  of  the  large  eddies  by  RANS,  gives  greater  generality,  in  principle, 
to  LES.  The  influence  of  the  unresolved  scales  of  motion  is  simulated  using  a  subgrid-scale  (SGS) 
model  (Smagorinsky  1963,  Lilly  1967,  Deardorff  1970,  Germano  et  al  1991,  Piomelli  et  al  1991, 
Ghosal  et  al  1995)  or  by  the  inherent  dissipation  in  the  numerical  scheme  (Boris  et  al  1992,  Oran 
and  Boris  1993,  Porter  et  al  1994,  Grinstein  1996,  Ansari  and  Strang  1996).  Because  the  statistics 
of  the  small  scale  turbulence  are  expected  to  be  more  homogeneous  and  isotropic  than  those  of  the 
large  scales,  a  general  model  of  the  small  scales  seems  more  plausible  than  a  general  model  of  the 
entire  spectrum  of  turbulent  motions. 

LES  has  been  shown  to  be  both  a  useful  research  tool  for  understanding  the  physics  of  turbulence, 
and  also  a  predictive  method  for  flows  of  engineering  interest.  Recent  compendia  and  reviews 
include  Galperin  and  Orszag  (1993),  Mason  (1994),  Lesieur  and  Metais  (1996)  and  Moin  (1997). 


Many  models  have  been  developed  for  the  subgrid-scale  stress  tensor.  These  include  the  conventional 
Smagorinsky  eddy  viscosity  model  (Smagorinsky  1963,  Lilly  1967,  Deardorff  1970),  the  spectra  eddy 
viscosity  model  of  Kraichnan  (1976),  the  dynamic  SGS  model  of  Germano  et  al  (1991),  the  scale 
similarity  model  of  Bardina  et  al  (1980),  and  the  localized  dynamic  SGS  model  of  Ghosal  et  al  (1995) 
and  more  recently  of  Menon  and  Kim  (1996),  and  many  others.  Although  most  research  has  focused 
on  incompressible  turbulent  flows,  there  has  recently  emerged  a  growing  interest  in  applications  of 
LES  to  compressible  turbulent  flows.  Examples  include  Yoshizawa  (1986),  Speziale  et  al  (1988), 
Moin  et  al  (1991),  Erlebacher  et  al  (1992),  Zang  et  al  (1992),  El-Hady  et  al  (1994),  Jansen  (1997), 
Spyropoulos  and  Blaisdell  (1996),  and  Haworth  and  Jansen  (1996).  Nearly  all  compressible  LES 
has  employed  spectral  methods  or  structured  grids,  with  the  exception  of  Jansen  and  Haworth. 

Apart  from  the  complexities  of  the  flowfield,  the  complicated  geometries  of  high  speed  vehicles  is 
also  a  challenge.  To  enable  treatment  of  complex  geometries  and  also  achieve  high  resolution  of 
the  flowfield  dynamically,  we  employ  an  unstructured  grid.  There  are  two  important  advantages 
of  unstructured  grids.  First,  algorithms  have  been  developed  to  facilitate  automatic  generation  of 
unstructured  grids  for  a  complex  geometries  (see,  for  example,  the  discussion  in  Barth  (1990,  1992). 
These  grid  generation  methods  can  be  substantially  more  efficient  (in  terms  of  user  time)  than  some 
of  the  multi-block  structured  grid  generation  methods  used.  Second,  local  mesh  refinement,  either 
adaptive  or  fixed,  can  been  performed  much  more  readily  for  unstructured  grids. 

The  report  summarizes  the  research  in  Large  Eddy  Simulation  of  compressible  turbulent  flows  using 
unstructured  grids.  Two  methods  for  simulation  of  the  subgrid  scale  stresses  have  been  examined. 
The  first  method  is  the  Monotone  Integrated  Large  Eddy  Simulation  (MILES)  technique.  The  sec¬ 
ond  method  is  a  hybrid  technique  combining  MILES  with  a  Smagorinsky  eddy  viscosity  model  for 
the  subgrid  scale  stresses.  These  two  methods,  together  with  different  algorithms  for  the  inviscid 
fluxes  and  function  reconstruction,  have  been  evaluated  for  four  turbulent  flows:  supersonic  bound¬ 
ary  layer,  supersonic  compression  corner,  supersonic  expansion-compression  corner  and  subsonic 
square  jet.  The  results  are  in  overall  good  agreement  with  the  experiment  and  Direct  Numerical 
Simulation  (DNS),  thereby  validating  the  accuracy  of  the  methodology. 

Governing  Equations 

The  governing  equations  are  the  three-dimensional  filtered  Navier-Stokes  equations.  For  a  function 
/,  its  filtered  form  /  is 

faHG,dr 

where  G  is  the  filtering  function,  and  its  Favre-averaged  form  /  is 


where  p  is  the  density.  From  the  Navier-Stokes  equations  for  the  instantaneous  flow  variables  density 
(p),  velocity  in  the  «th  coordinate  direction  (itj),  pressure  (p)  and  temperature  (T),  Favre-averaging 
and  spatial  filtering  yield  the  filtered  Navier-Stokes  equations  (here  written  using  the  Einstein 
summation  notation  where  repeated  indices  denote  summation) 
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Two  different  Sub-Grid-Scale  (SGS)  models  are  employed.  The  first  model  is  Monotone  Integrated 
Large  Eddy  Simulation  (MILES)  wherein  the  numerical  algorithm  itself  provides  the  requisite 
dissipation  associated  with  the  subgrid  scale  motions.  The  second  model  is  the  classical  constant- 
coefficient  Smagorinsky  method 


Sij 


Tij 

Qj 


1  ( duj  diij\ 

2  \dxj  dxi  J 

2CrpA  \J Smn Smn  (§ij  3  Skk^ij 


Pcp 


Pr, 


df 

dxj 


where  Cr  =  0.00423  and  A  is  the  length  scale  which  is  related  to  the  local  grid  size.  For  boundary 
layer  flows,  A  is  multiplied  by  the  Van  Driest  damping  factor 

D  =  1  -  e~n+'A 

where  A  —  26,  n+  =  nuT/vw  is  the  normal  distance  to  the  (nearest)  solid  boundary  normalized  by 
the  viscous  length  scale  vw/uT  where  vw  is  the  kinematic  viscosity  evaluated  at  the  wall  and  UT  is 
the  local  friction  velocity. 

We  simplify  the  notation  by  hereafter  dropping  the  tilde  ~  and  over  bar  ~.  The  flow  variables  are 
nondimensionalized  using  the  reference  density  Pco,  velocity  Uoo,  static  temperature  and  length 
scale  L,  with  Mach  number  =  Uco/y/,yRTrx>.  The  governing  equations  are  therefore 
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Numerical  Algorithm 

The  governing  equations  are  expressed  in  finite  volume  form  for  a  control  volume  V  with  surface 
dV 

4 :[QdV+[  {Ft  +  Gj  +  Hk)  •  ndA  —  0 
at  Jv  JdV 


where  Q  is  the  vector  of  dependent  variables 


Q  = 


P  \ 
pu 
pv 
pw 

\  Pe  / 


and  the  flux  vectors  are 


pu  \ 

pv  \ 

f  pw  \ 

pu 2  +  p  -  Txx 

puv  -  TXy 

puw  -  Txz 

puv  -  Txy 

,  G  = 

PV 2  +  P  —  lyy 

,  H  = 
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{pe  +  p)u  -  Qx  -  px  ) 

{pe+p)v  -Qy-Py  ^ 

\  (pe  +  p)w  -  Qz-  Pz  / 

with 

fix  =  '~J~xx'Uj  'Txy'V  “f“  'Txz'W 
fiy  =  TXyU  -f*  'TyyV  ~h  'Tyz'W 
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An  unstructured  grid  of  tetrahedra  is  employed,  with  a  cell-centered  storage  architecture.  The 
cell-averaged  values,  stored  at  the  centroid  of  each  tetrahedron  of  volume  V\  are 

Qi  =  ^  [  Qdv 

Vi  JVi 

The  inviscid  fluxes  are  computed  using  Godunov’s  method  which  is  an  exact  one-dimensional  Rie- 
mann  solver  (Gottlieb  and  Groth  1988)  applied  normal  to  each  face.  The  inviscid  flux  computations 
®  require  the  values  of  each  variable  on  either  side  of  the  cell  faces.  These  values  are  obtained  from  the 

cell-averaged  values  by  second-order  or  third-order  function  reconstruction  using  the  Least  Squares 
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method  of  Ollivier-Gooch  (Ollivier-Gooch  1997).  The  second-order  function  reconstruction  method 
of  Prink  (1994)  was  employed  in  some  of  the  earlier  LES  studies,  but  was  found  inferior  to  the 
method  of  Ollivier-Gooch  (Okong’o  and  Knight  1998).  More  details  on  the  reconstruction  schemes 
are  given  in  Okong’o  and  Knight  (1998). 

The  viscous  fluxes  and  heat  transfer  are  computed  by  application  of  Gauss’  theorem  to  the  control 
volume  whose  vertices  are  the  centroids  of  the  cells  which  share  each  node.  The  second-order 
accurate  scheme  (in  2-D)  is  given  by  Knight  (1994)  and  the  extension  to  3-D  is  straightforward. 

Parallelization 


The  code  is  parallelized  using  domain  decom¬ 
position  and  Message  Passing  Interface  (MPI). 
Domain  decomposition  is  performed  in  a  pre¬ 
processing  step.  The  domain  is  decomposed  in 
a  single  direction  with  equal  number  of  tetrahe- 
dra  in  each  domain.  A  halo  of  cells  is  added  in 
each  domain  to  provide  data  on  the  adjacent  do¬ 
main,  and  the  halo  cell  data  is  updated  at  ev¬ 
ery  subiterate  of  the  time  integration.  An  ex¬ 
ample  is  shown  in  Fig.  1  for  the  LES  of  decay 
of  isotropic  turbulence.  The  numerical  algorithm 
achieves  excellent  parallel  performance.  For  ex¬ 
ample,  the  speed-up  on  four  processors  of  the 
SGI  Power  Onyx  with  R-10000  processors  is  3.7 
for  93%  efficiency  (Knight  et  al  1998). 


Figure  1:  Example  of  domain  decomposition 


Results 

Four  different  configurations  have  been  examined:  supersonic  flat  plate  boundary  layer,  supersonic 
compression  corner,  expansion-compression  corner  and  subsonic  square  jet. 


1  Supersonic  Flat  Plate  Turbulent  Boundary  Layer 

The  adiabatic  and  isothermal  flat  plate  turbulent  boundary  layers  at  Mach  3  and  Mach  4  at 
Reynolds  number  Reg  =  2  x  104  (based  on  the  incoming  boundary  layer  thickness  5)  have  been 
computed.  The  Reynolds  number  based  on  the  momentum  thickness  62  and  wall  viscosity  yw  is 
Re, 52  =  600.  The  Reynolds  number  is  sufficently  high  to  achieve  turbulent  flow. 

The  inflow  conditions  are  obtained  using  a  compressible  extension  of  the  method  of  Lund  et  al 
(1998).  The  simulation  generates  its  own  inflow  conditions  through  a  sequence  of  operations  where 
the  velocity  field  at  a  downstream  station  is  rescaled  and  reintroduced  at  the  inflow  boundary 
(Fig.  2).  Defining  x,y  and  z  to  denote  the  streamwise,  transverse  and  spanwise  directions,  re¬ 
spectively,  the  size  of  the  computational  domain  is  Lx  =  14.85,  Ly  =  3.45  and  Lz  =  2.05.  The 
spanwise  width  Lz  is  approximately  three  times  the  experimental  spanwise  streak  spacing  (as¬ 
suming  the  compressible  turbulent  boundary  layer  streaks  scale  in  accordance  with  incompressible 
experimental  results) .  The  streamwise  length  Lx  is  approximately  three  times  the  mean  experimen- 


Figure  2:  Computational  domain 


tal  streamwise  streak  size.  The  height  Ly  is  based  on  the  requirement  that  acoustic  disturbances 
originating  at  the  upper  boundary  do  not  interact  with  the  boundary  layer  on  the  lower  wall. 

The  reference  quantities  for  non-dimensionalization  are  the  incoming  boundary  layer  thickness  £, 
velocity  [Too,  density  static  temperature  Too  and  molecular  viscosity  p0 0  (where  the  subscript  oo 
denotes  the  freestream  condition).  The  grid  resolution  near  the  wall  is  dependent  on  Ax+,  A y+  and 
A 2+,  where  Ax+  =  Ax/77,  Ay+  “  A77/77  and  A z+  =  Az/rj.  The  inner  length  scale  is  77  =  vw/uT, 
where  uw  is  the  kinematic  viscosity  at  the  wall,  uT  =  \frw/pw  is  the  friction  velocity,  rw  is  the  wall 
shear  stress  and  pw  is  the  density  at  the  wall.  Before  we  proceed  with  the  discussion  about  the  grid 
resolution,  we  first  describe  how  to  obtain  77. 

The  theoretical  value  of  the  friction  velocity  uT  for  a  supersonic  flat  plate  boundary  is  obtained 
from  the  combined  Law  of  the  Wall  and  Wake  evaluated  at  y  =  5 


Uvd  l  ,  ur  2U  .  2 

- =  —  ln(y — )  +  C  + — sin  ““ 

Ur  K  Vn)  K  2  0 


(l) 


where 


U 


UVD  =  {sin-1! 


ir2A2(^)-5 


VB2  +  4A2 


VB2  +  4A2 


]} 


A 

B 


=  J'oo  (^T“)  + 

1 00 

=  (^PUrnMl^2 

__  Taw 
T^~ 

=  Too[l  + 


(2) 


where  k  =  0.4  is  von  Karman’s  constant,  C  —  5.1,  the  wake  parameter  n  is  0.12  at  Reg  =  2  x  104, 
the  exponent  u  is  0.76,  the  mean  turbulent  Prandtl  number  Prtm,  is  0.89  and  the  ratio  of  specific 
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heats  7  is  1.4.  The  wall  temperature  Tw  is  fixed  at  10%  above  the  theoretical  adiabatic  temperature 
Taw  for  the  isothermal  boundary.  In  the  computation,  uT  and  uw  are  obtained  from  uT  =  \Jtv,/ pw 
and  vw  =  ^oo(y^)1+a;,  respectively. 

The  variation  of  Mach  number  and  the  different  temperature  boundary  condition  have  effect  on 
7,  therefore  on  Ax+,  Ay+  and  A z+.  At  the  same  Mach  number,  the  wall  temperature  for  the 
isothermal  case  is  10%  higher  than  that  for  the  adiabatic  case,  leading  to  the  larger  y  and  the 
smaller  Ax+  if  keeping  the  same  Ax.  However  this  effect  is  very  small  in  our  case,  therefore  we 
keep  the  same  Ax,  Ay  and  Az  for  the  different  temperature  condition  at  the  same  Mach  number. 
The  grid  details  are  shown  in  Table  2. 


Table  2:  Details  of  Grid 


A3 

13 

A4 

14 

Ax+ 

20 

18 

12 

11 

A  y+ 

1.8 

1.6 

1.8 

1.6 

A  z+ 

7 

6.4 

4 

3.4 

Ax/8 

0.1 

0.1 

0.1 

0.1 

Ay/8 

0.18 

0.18 

0.18 

0.18 

Az/8 

0.034 

0.034 

0.034 

0.034 

where  A  and  I  stand  for  the  adiabatic  and  isothermal  cases,  respectively  and  the  number  followed 
indicates  the  Mach  number.  The  Ax+,A y+  and  A z+  are  measured  at  the  wall  and  the  Ax,  Ay 
and  Az  are  measured  at  y/S  =  1.0.  The  grid  is  uniform  in  x  and  z  directions  and  stretched  in  y 
direction  with  about  23  layers  of  tetrahedra  in  the  boundary  layer  for  each  case. 

The  initial  condition  is  a  turbulent  mean  profile  with  random  fluctuations.  The  simulation  is  run 
first  for  90  inertial  timescales  S/Uoo  in  order  to  eliminate  starting  transients  (Lund  et  al  1998). 

For  a  function  /,  its  average  in  time  form  <  /  >  is  defined  by 

1  rlf 

</>=  7 - t  fdt 

tf  -  H  Jti 

and  its  time  fluctuating  part  is 

/"  =  /-</> 


In  order  to  provide  converged  data,  the  primitive  variables  are  averaged  in  spanwise  direction  and 
the  statistical  evaluations  are  performed  on  a  period  longer  than  tf  —  t{  =  4(W/t/oo-  The  notation 
for  the  combined  temporal  and  spanwise  average  is 

1  1  fLz  flf 

</»=———/  /  fdtdz 

L/z  tf  H  J 0  Jti 

A  simplifying  notation  is  used  for  the  velocity,  temperature  and  pressure 

U  =<C  u  » 


The  mean  streamwise  velocity  profiles  using  the  Van  Driest  transformation  are  plotted  in  Fig.  3 
and  Fig.  4  (where  uT  is  obtained  from  the  simulation).  Good  agreement  is  shown  with  the  viscous 
sublayer  linear  approximation  Uvd/ut  =  y+  and  Law  of  the  Wall  formulated  in  (1). 


Figure  3:  Van  Driest  velocity  at  M=3 


Figure  4:  Van  Driest  velocity  at  M=4 


The  mean  velocity  profiles  shown  in  Fig.  5  and  Fig.  6  exhibit  virtually  identical  distributions 
for  adiabatic  and  isothermal  cases  and  show  good  agreement  with  experiment  (Zheltovodov  et 
al  1986,  Zheltovodov  et  al  1990).  The  mean  temperature  profiles  in  Fig.  7  and  Fig.  8  display  a 
higher  temperature  distribution  for  the  isothermal  case  with  the  wall  temperature  higher  than  the 
adiabatic  and  experiment  data  (Zheltovodov  et  al  1990)  which  are  also  obtained  at  the  adiabatic 
boundary  condition.  The  difference  is  expected  since  the  wall  temperature  for  the  isothermal  case 
is  fixed  at  10%  higher  than  the  adiabatic  case. 


Figure  5:  Streamwise  velocity  at  M=3  Figure  6:  Streamwise  velocity  at  M=4 

The  discrepancies  between  Mach  4  cases  and  experiments  in  Fig.  6  and  Fig.  8  are  due  to  the  effects 
of  Mach  number  and  Reynolds  number.  The  Reynolds  number  in  the  simulation  is  one  magnitude 
lower  than  experiments  due  to  the  significant  computation  cost  in  LES.  The  outer  portion  of  the 
velocity  profiles  in  Fig.  6  is  in  good  agreement  with  experiment  since  this  portion  is  not  sensitive 
to  the  Reynolds  number.  The  discrepancy  in  the  inner  portion  is  due  to  the  effect  of  the  Reynolds 


Figure  7:  Temperature  at  M~ 3 


Figure  8:  Temperature  at  M=4 


number.  The  effect  of  Mach  number  is  observed  in  Fig.  8,  which  can  be  explained  using  Crocco’s 
relationship  between  the  mean  temperature  and  mean  velocity  profiles 


T  Tw  Taw  Tw  U  T  ^  ^  \2 

=  j*  *"  — r — TT  ”  r”2“Moo( 

-*00  1 OO  -*-00  uoo  6  u  oo 

where  r  is  the  recovery  factor  defined  as 

__  Tr  ~  T0 Q 
"'“To-Too 


(3) 

(4) 


where  Tr  is  the  adiabatic  or  recovery  temperature  and  To  is  the  freestream  stagnation  temperature. 
For  the  adiabatic  case,  Eq.  (3)  becomes 


OO 


(5) 


Eqs.  (3)  and  (5)  show  the  trend  that  under  the  same  mean  velocity  distribution,  the  mean  temper¬ 
ature  decreases  with  increasing  Mach  number,  leading  to  the  discrepancy  between  the  calculation 
and  experiment  in  the  outer  portion  of  the  mean  temperature  profiles  in  Fig.  8.  The  discrepancy 
in  the  inner  portion  is  mainly  due  to  the  effect  of  the  Reynolds  number. 

The  mean  streamwise  resolved  turbulent  kinematic  normal  stress  C  unun  »,  normalized  using  the 
local  mean  density  and  wall  shear  stress  rW:  is  shown  in  Fig.  9  and  Fig.  10.  As  discussed  in 

Zheltovodov  and  Yakovlev  (1986)  and  Smits  and  Dussauge  (1996),  the  scaling  <§C  p  unu11  > 
/ rw  provides  an  approximate  self-similar  correlation  of  experimental  data  for  supersonic  flat  plate 
zero  pressure  gradient  adiabatic  boundary  layers,  although  the  measurements  close  to  the  wall  are 
subject  to  considerable  uncertainty.  In  those  two  figures  data  are  displayed  from  Konrad  and  Smits 
(1998),  Johnson  and  Rose  (1975),  Muck  et  al  (1984,  1985),  Konrad  (1993),  as  well  as  upper  and 
lower  bounds  of  an  extensive  set  of  experimental  data  for  the  Mach  number  range  M  =  1.72  to 
9.4  in  accordance  with  generalizations  of  Zheltovodov  and  Yakovlev  (1986).  The  characteristics  of 
the  different  experiments  are  displayed  in  Table  3.  The  computed  results  show  good  agreement 
with  experiment  for  the  main  part  of  the  boundary  layer  ( y/S  >  0.2),  despite  a  significantly  higher 


Table  3:  Flat  Plate  Boundary  Layer  Experimental  Data 


Name 

Mach  No. 

Res 

LES 

3.0  &  4.0 

20  x  103 

DNS  Adams  (1997) 

3.0 

25  x  103 

Johnson  &  Rose  (1975) 

2.9 

1000  x  103 

Konrad  (1993) 

2.9 

1590  x  103 

Konrad  &  Smits  (1998) 

2.87 

1900  x  103 

Muck  et  al  (1984,  1985) 

2.87 

1638  x  103 

Zheltovodov  et  al  (1986) 

1. 7-9.4 

up  to  2000  x  103 

experimental  Reynolds  number.  The  decreasing  slope  corresponds  precisely  to  Johnson  and  Rose 
(1975)  data.  For  y/S  <  0.2  the  presence  of  the  typical  high  level  peak  in  the  near  wall  region  is 
#  supported  by  experimental  data  of  Konrad  (1993)  and  the  Direct  Numerical  Simulation  data  from 

Adams  (1997),  which  is  nearly  at  the  same  Reynolds  number  as  the  LES.  However,  no  conclusion 
can  be  drawn  about  the  precise  y  position  and  the  width  of  this  peak  without  further  experimental 
data  or  DNS. 


Figure  9:  Streamwise  Reynolds  stress  at  M=3  Figure  10:  Streamwise  Reynolds  stress  at  M=4 

In  Fig.  11  and  Fig.  12  Reynolds  shear  stress  distributions  are  shown  for  the  same  experiments  and 
the  DNS.  Again,  the  data  fit  well  in  the  outer  part  of  the  boundary  layer.  The  maximum  value  and 
the  decreasing  slope  are  again  well  predicted. 


The  capability  of  our  LES  method  to  accurately  predict  the  heat  transfer  in  the  flat  plate  boundary 
layer  is  evaluated.  The  Reynolds  analogy  relates  the  skin  friction  coefficient  Cf  and  heat  transfer 
coefficient  Ch  by  the  Prandtl  number  as  follows 


where  Ch  and  Cf  are  written  as 


2  Ch  1 
Cf  Pt  tin 


Ch  = 


Qw 


P<x>Uoc>Cp{Tw  Taw) 


(6) 

(7) 
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Figure  11:  Reynolds  shear  stress  at  M=3 


Figure  12:  Reynolds  shear  stress  at  M=4 


;  boo  ui 

where  cp  is  the  specific  heat  at  constant  pressure.  The  wall  heat  flux  (qw)  and  skin  friction  (tw)  are 
obtained  from  the  isothermal  case  and  the  adiabatic  wall  temperature  ( Taw )  is  calculated  from  the 
adiabatic  case.  The  wall  heat  flux  is 

qw  =  —A  (9) 

dy 


and  the  wall  shear  stress  is 


Tw  —  " 


Table  4:  LES  predictions 


Name  T^/Tx.  C/  Ch 

A2.88  2^51  2.44  x  10~3  0 

12.88  2.72  2.32  x  10“3  1.27  x  1CT3 

~M  T95  1.97  x  ltH  0 


2.00  x  10- 


1.24  x  10“ 
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Table  5:  Comparison  of  LES  and  Experiment 


Cases 

Name 

LES 

Experiment 

Error 

Mach=2.88 

Taw /Too 

2.51 

2.549 

1.5% 

Cf 

2.32  x  1(T3 

2.56  x  lO"3 

9.4% 

Ch 

1.27  x  Hr3 

1.44  x  lO"3 

11.8% 

2Ch/Cf 

1/0.91 

1/0.89 

2.2% 

P^tm 

0.91 

0.89 

2.2% 

Mach=4.0 

Taw  /Too 

3.95 

3.848 

2.6% 

Cf 

2.0  x  lO-3 

2.17  x  lO-3 

7.8% 

Ch 

1.24  x  lO-3 

1.22  x  lO-3 

1.6% 

2Ch/Cf 

1/0.81 

1/0.89 

9.9% 

P^tm 

0.81 

0.89 

9.0% 

Note: 


Experimental  Taw  from  Eq.  (3) 
Experimental  Cf  from  Eqs.  (1),  (2)  and  (8) 
Experimental  Ch  from  Eq.  (6) 


The  wall  temperature  Tw  is  fixed  at  Tw  =  1 .  lTaw ,  where  the  empirical  adiabatic  wall  temperature 
Taw  =  [1  +  PrtmM%o]Too.  All  the  predicted  results  are  listed  in  Table  4,  where  the  wall 
temperature  for  the  isothermal  case  is  fixed.  The  comparison  with  experiment  is  shown  in  Table  5. 
The  experimental  adiabatic  wall  temperature  is  computed  from  Taw  =  [1  +  ^Y^PrtmM^Too.  The 
computed  mean  turbulent  Prandtl  number  from  (6)  shows  good  agreement  with  experiment  value 
of  0.89,  indicating  the  consistency  of  LES  results  with  the  Reynolds  analogy. 


The  turbulent  Prandtl  number  changes  across  the  boundary  layer.  Simpson  et  al  (1970)  have  estab¬ 
lished  the  uncertainty  envelope  of  the  turbulent  Prandtl  number  for  incompressible  zero  pressure 
gradient  turbulent  boundaries.  The  experimental  predictions  by  Meier  and  Rotta  (1971)  at  Mach 
number  up  to  4.5  at  the  wall  and  Horstman  and  Owen  (1972)  at  M=7.2  and  cooled  wall  conditions 
fall  into  this  uncertainty  envelope.  According  to  the  eddy  viscosity  hypothesis,  the  turbulent  stress 
and  heat  flux  can  be  expressed  as 


dui  duk 


Tik  =  -pu'jU'f.  =  Mt(^—  + 


dx k  dxi 


'idik )  -  -pkdik 

3  dxj  ’  y 


(11) 


Qk  —  cppT'u'k  —  A  (12) 

where  the  bar  “denotes  the  filtered  flow  variables  and  the  tilde 'denotes  the  Favre- averaged  filtered 
flow  variables.  The  local  turbulent  Prandtl  number  ( Prt )  for  a  two-dimensional  boundary  layer  can 
be  derived  from  the  above  two  equations  as 


Pr 


_  % 
t  du 
dy 


pT'v' 


(13) 


The  calculated  turbulent  Prandtl  number  profile  is  shown  compared  with  the  experimental  range 
in  Fig.  13.  The  Prandtl  number  reaches  the  maximum  at  the  wall  and  starts  to  decrease  away  from 
the  wall.  In  the  outer  portion  of  boundary  layer,  the  fluctuation  of  Prandtl  number  is  relatively 
greater  than  the  inner  portion,  which  is  consistent  with  the  experimental  trend. 
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Figure  14:  Example  of  overshoot 


Figure  13:  Turbulent  Prandtl  number 

2  Compression  Corner 

Supersonic  flow  past  a  compression  corner  is  an  important  problem  in  aerodynamics.  It  represents, 
for  example,  the  deflection  of  a  control  surface  on  a  wing.  The  shock  can  cause  a  boundary  layer 
separation  upstream  of  the  point  of  impingement  of  the  primary  shock,  with  a  secondary  shock 
forming  near  the  separation  (Andreopoulos  and  Muck  1987,  Dolling  and  Or  1983,  Horstman  et  al 
1977,  Settles  et  al  1979,  Smits  and  Muck  1987,  Zheltovodov  et  al  1983,  Zheltovodov  and  Yakolev 
1986,  and  Zheltovodov  1996).  Reynolds-averaged  Navier-Stokes  simulations  have  failed  to  accurately 
predict  the  flow  characteristics  (Knight  and  Degrez  1998)  such  as  fluctuating  pressure  and  heat 
transfer. 

2.1  Weighting  Function  and  Limiters 

The  computation  of  a  strong  shock  using  the  exact  Riemann  solver  (Godunov’s  method)  and  the 
Least  Squares  method  (Ollivier  1997)  as  a  reconstruction  scheme  leads  to  the  generation  of  oscilla¬ 
tions  in  the  vicinity  of  shock  waves  which  cause  numerical  instability.  It  can  be  shown  theoretically 
that  a  linear  second-order  upwind  scheme  always  generates  oscillations.  The  only  way  to  overcome 
this  limitation,  while  satisfying  the  concept  of  an  entropy  function,  is  to  introduce  non-linear  com¬ 
ponents.  The  classical  method  to  avoid  such  spurious  oscillations  is  to  implement  limiters  which 
control  the  gradient  of  the  computed  quantities  to  prevent  the  appearance  of  overshoots  and  under¬ 
shoots.  An  excellent  review  of  this  technique  is  described  in  Hirsh  1997.  A  different  approach  is  the 
stencils-searching  ENO  schemes  which  have  been  extended  to  unstructured  grids  (Abgrall  1994). 
At  each  time  step,  the  stencil  is  chosen  which  provides  the  smoothest  reconstruction.  However, 
this  approach  is  computationally  expensive  for  LES  since  it  implies  a  determination  of  the  stencil 
at  every  time  step.  Recently,  Ollivier-Gooch  (Ollivier  1997)  proposed  a  weighted  stencil  method 
wherein  the  stencil  is  fixed  but  the  weights  are  recomputed  at  each  time  step  as  required1. 

xWe  found  that  the  weighted  stencil  method  of  Ollivier-Gooch  (Ollivier  1997)  provided  an  improvement  compared 
to  the  unweighted  results,  but  was  nonetheless  very  sensitive  to  the  weighting  parameters.  For  the  25°  compression 
corner,  we  found  that  overshoots  and  undershoots  in  the  vicinity  of  the  shock  could  not  be  avoided  without  adversely 


2.2  Limiters 


We  consider  limiters  which  control  the  gradient  of  computed  quantities  reconstructed  to  a  cell  face 
(denoted  by  the  index  i  +  5).  The  limiters  described  in  the  literature  ( e.g Van  Leer’s  limiter, 
Minmod,  Roe’s  Superbee  limiter)  are  expressed  as  a  function  of  the  ratio  ri+ 1  of  consecutive 
variations  2 


r- 


i+i 


Uj+l  Uj 
Ui  Ui—i 


(14) 


This  expression  is  well  defined  in  case  of  a  structured  grid,  where  i  +  1  means  the  next  cell  in  the 
i  discretization  and  i  —  1  means  the  previous  one. 


2.3  Homogeneous  Limiter 


Consider  a  linear  reconstruction  of  a  flow  variable  wherein  the  computed  gradient  yields  an  over¬ 
shoot  (or  undershoot)  in  the  reconstructed  value  at  the  cell  interface  (Fig.  14).  The  overshoot  (or 
undershoot)  can  be  avoided  if  the  interface  values  were  to  remain  between  the  adjacent  cell  aver¬ 
aged  values.  We  can  limit  the  slope  computed  by  the  LS  method  to  satisfy  this  criterion.  For  a  one 
dimensional  case,  the  interface  value  computed  from  the  left  cell  i  is 


=Ui  +  Ci(xi+ 1  -  Xi) 


(15) 


where  C{  is  du/dx.  Consider  the  case  Ui  <  ui+1.  Referring  to  Fig.  14,  the  reconstructed  value  ulei\ 

M "  2 

is  required  to  lie  within  the  adjacent  cell  averaged  values 


Ui  <  u*!\  <  ui+ 1 


or 


0  <Ci< 


'U’i 

xi+i  —  xi 

This  is  achieved  by  replacing  the  gradient  C{  by  TjCi  where 

1 


— ni-  \ 


for  0  <  Ct  < 

-  "i+j-*.- 
r<  ^  Ui+i-Ui 
for  c' 


e‘+r*1 


0 


for  Ci  <  0 


An  analogous  result  can  be  obtained  for  the  case  u ;  >  u  • ,  i . 

l'  2 

For  a  general  3D  configuration,  the  limited  quantities  are 


left  _ 


Ui,3 


=  Ui  +  r)C  ■  Aa; 


(16) 

(17) 


(18) 


(19) 


where  ufj*  is  the  reconstructed  value  of  variable  u  for  cell  i  and  the  face  adjacent  to  cell  j,  and 


C- Ax  = 


Cxi  \ 

Cyi 

Czt  / 


^  xi,j 
Vi, 3  ~  Vi 
\  zi,j  ~  zi 


,  \ 
/ 


(20) 


affecting  the  undisturbed  boundary  layer. 


Figure  15:  Temperature  (no  limiter) 


Figure  16:  Temperature  (homogeneous  limiter) 


V- 


1.0 

for 

0  <  C  •  Ax  <  Uj  —  U{ 

Uj  —Ui 

for 

C  •  Ax  >  Uj  ~  Uj 

C- Ax 

J  b 

— #  — * 

0 

for 

C  ■  Az  <  0 

(21) 


The  above  expression  holds  for  U{  <  Uj.  An  analogous  expression  holds  for  Uj  >  uj.  In  practice  the 
limiter  is  successively  applied  to  the  different  cell  neighbors  j  of  the  cell  i .  As  far  as  tetrahedras  are 
concerned,  only  the  neighbors  sharing  a  face  are  used. 


A  25°  compression  corner  computation  has  been  designed  to  evaluate  the  efficiency  of  the  limiter. 
Fig.  15  displays  the  evolution  of  the  static  temperature  across  the  shock  computed  using  the  second 
order  LS  without  a  limiter.  Strong  overshoots  and  undershoots  appear.  Fig.  16  displays  the  static 
temperature  in  the  shock  using  this  limiter.  The  spurious  oscillations  disappear,  although  the 
gradients  within  the  cells  appear  to  have  been  reduced  probably  more  than  necessary. 


2,4  Inhomogeneous  Limiter 


The  homogeneous  limiter  (21)  does  not  treat  the  function  gradients  Cxj,Cyi,Czi  independently, 
and  therefore  is  overly  limiting.  Consider  a  case  where  Cx{  and  Cyi  are  “satisfactorily”  computed 
but  Czi  is  overestimated  by  the  LS  method  and  therefore  ufj1  is  also  overestimated.  During  the 

limitation  process  rj  is  set  to  a  value  less  than  1  because  C  •  Ax  is  too  large.  As  a  result,  all  three 
gradients  are  reduced,  even  though  Cx,Cy  were  “satisfactory”. 

A  three  parameter  limiter  is  defined  by 


u, 


left  _ 


h3 


—  U{  + 


1 lxCxi  ^ 

r]yCyi 

vzCzi  y 


X{j  X{ 

Uij  ~~  Vi 

ZiJ  -  Zi 


(22) 


wherein  each  gradient  Cxi,  Cyi ,  Cz{  has  a  limiter.  We  present  the  inhomogeneous  limiter  for  Ui  <  Uj. 
Analogous  expressions  hold  for  Uj  >  uj. 
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— *  — * 

Case  1:  C  •  Ax  >  Uj  —  ui 


tto  =  1.0-fc1|C^-^| 
%  =  i.o-fci|Qs^pal| 


where 


*1 


~left 

_  Uu  ~ 

*left 

uJ  ~ui 


Cxi.(xij  -  Xi ) 
k2  =  max  <  Cyi.(yij  -  yt) 

.  Czi.(zij-Zi) 

where  is  the  reconstructed  value  assuming  T)x  =  r)y  =  )jz  =  1 


Cose  0  <  C  •  Aa:  <  Uj  —  U{ 


(23) 


(24) 

(25) 


Case  5:  C  •  Ax  <  0 


%  =  1 
%  =  1 
=  1 


r)x  =  o 

=  o 

^2  =  0 


(26) 

(27) 

(28) 


(29) 

(30) 

(31) 


The  process  mainly  limits  the  components  rix,riy,r)z  whose  contribution  in  the  C  ■  Ax  expression  is 
the  largest.  When  the  j  neighbor  position  is  aligned  with  x,y  ov  z  direction,  then  only  Cx,  Cy  or 
Cz  is  respectively  reduced. 

Fig.  17  displays  the  results  of  this  limiter  for  the  static  temperature  profile  within  the  shock.  The 
effect  on  the  gradient  is  less  severe  than  in  Fig.  16.  The  25°  compression  corner  computations 
described  below  use  this  limiter. 


^  2.5  ENO  scheme 

The  key  idea  of  ENO  schemes  is  to  use  the  smoothest  stencil  among  several  candidates  to  approxi¬ 
mate  the  fluxes  at  cell  boundaries  to  a  high  order  accuracy  and  at  the  same  time  to  avoid  spurious 
oscillations  near  shocks.  We  previously  computed  supersonic  turbulent  flow  past  compression  cor¬ 
ner  of  8°  at  Mach  3  without  use  of  limiters  (Urbin  1999).  However,  the  25°  compression  corner  at 
®  Mach  3  (see  Fig.  18)  requires  a  limiter,  and  thus  determination  of  the  most  effective  limiter  is  the 

objective  of  this  part  of  our  study.  In  the  25°  compression  corner,  a  strong  shock  is  formed  through 
a  series  of  unsteady  compression  waves,  and  variations  in  the  flow  variables  in  the  boundary  layer 
are  sometimes  comparable  to  the  shock.  In  order  to  determine  a  proper  criterion  to  construct  the 
ENO  stencil,  we  begin  with  the  analysis  of  the  variable  gradient  (we  use  density  since  it  has  a 
#  similar  change  as  the  other  variables),  shown  in  Fig.  19  at  x  =  3 S,  where  x  is  measured  from  the 

corner  and  6  is  the  inflow  boundary  layer  thickness.  The  density  is  non-dimensionalized  by  the 
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Figure  17:  Temperature  (inhomogeneous  limiter)  Figure  18:  25°  compression  corner 


freestream  density  p^.  The  gradient  using  the  isotropic  stencil  (Okongo  1998)  with  Least  Squares 

9  reconstruction  shows  two  peaks  corresponding  to  the  boundary  layer  and  the  shock  wave.  Outside 

these  two  regions,  the  change  of  the  gradient  is  smooth.  The  selection  of  the  ENO  stencil  should 
take  advantage  of  this  feature  and  make  the  ENO  stencil  focus  on  these  two  regions. 

First,  we  establish  a  criterion  to  decide  where  to  use  the  ENO  stencil.  The  density  ratio  across  the 
shock  is  approximately  2.19  according  to  the  Rankine-Hugoniot  equations.  The  maximum  norm  of 

#  the  density  gradient  a  is  approximately  equal  to  6  with  mesh  spacing  A r  =  \J Aa:2  +  Ay2  +  Az2  = 
0.2  in  the  vicinity  of  the  shock,  where  A r  is  non-dimensionalized  by  S.  We  make  a  cut-off  at  a  =  6 
and  any  cell  whose  density  gradient  is  larger  than  a  should  use  an  ENO  stencil.  We  define  E  as  a 
set  of  cells  satisfying  |Vp|  >  a.  In  Fig.  19,  the  maximum  density  gradient  based  on  the  isotropic 
stencil  is  approximately  13. 

♦  Second,  for  all  E,  a  direction  should  be  found  to  construct  the  ENO  stencil.  Consider  any  tetrahe¬ 
dron  in  E  in  which  four  possible  density  gradients  are  computed  using  only  three  face  neighbors, 
denoted  as  |Vp|*  ( k  =  0,1,2, 3).  The  subscript  is  the  index  of  the  face  neighbor  which  has  been 
excluded.  The  density  gradients  are  computed  by  the  Least  Squares  method.  We  assume  that 


|Vp|3  >  |Vp|2  >  |Vp|i  >  |Vp|o 
where  |Vp|max  =  |Vp|3,  |Vp|min  =  |Vp|0. 

We  construct  the  ENO  stencil  as  follows.  First,  the  three  face  neighbors  which  give  the  minimum 
density  gradient  are  added.  Second,  we  find  all  the  node  neighbors  of  any  cell  of  E  and  exclude  all 
the  cells  which  share  at  least  one  node  with  the  remaining  face  neighbor.  Any  five  of  the  remaining 
cells  which  are  the  node  neighbors  of  that  cell  and  these  three  face  neighbors  are  used  to  construct 
the  ENO  stencil. 


A  modified  Riemann  shock  tube  case  with  the  initial  pressure  ratio  of  ten  and  an  initial  distribution 
of  isotropic  turbulence  is  used  to  verify  the  ENO  scheme  and  to  compare  with  the  inhomogeneous 
limiter.  The  pressure  ratio  across  the  shock  is  chosen  close  to  that  in  the  25°  compression  corner. 
The  cut-off  number  a  is  set  to  6.  Fig.  20  shows  the  comparison  between  the  inhomogeneous  limiter 
and  the  ENO  scheme  for  the  Riemann  shock  tube  case.  It  can  be  seen  that  the  ENO  stencil  is 
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Figure  19:  |Vp|  vs  y/S 
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Figure  21:  ENO  cells  for  25°  compression  corner  Figure  22:  Flowfield  behaviour  across  shock 


used  only  in  the  shock  region,  while  the  inhomogeneous  limiter  is  also  used  in  the  expansion  fan 
region  where  there  are  no  overshoots  and  even  in  the  undisturbed  flow  in  presence  of  turbulence. 
In  Fig.  21,  the  ENO  scheme  is  applied  to  the  25°  compression  corner  with  a  =  6.  The  cells  using 
the  ENO  stencil  appear  principally  in  the  shock  region  and  the  boundary  layer  downstream  of 
the  corner  where  the  density  gradient  has  a  larger  variation  than  the  upstream.  Fig.  22  shows  the 
instantanous  profiles  of  velocity  and  pressure  along  the  direction  perpendicular  to  the  shock  wave 
at  z  =  1.0<5  (where  ye  is  measured  along  the  direction  perpendicular  to  the  shock  wave  and  Uoo  is 
the  streamwise  velocity  in  the  freestream.).  The  x  coordinate  of  the  interception  point  of  the  ye 
with  the  solid  wall  is  5<5,  where  x  is  measured  from  the  corner.  The  shock  wave  is  efficiently  limited 
within  the  width  of  two  or  three  cells. 

2.6  Results 

David  (1993)  performed  the  first  LES  of  a  compression  corner  which  successfully  reproduced  the 
Taylor-Gortler  vortices  downstream  of  the  shock.  Nevertheless,  the  use  of  a  pseudo-compressible 
subgrid-model  did  not  permit  accurate  quantitative  results.  The  second  and  most  recent  LES  was 
Hunt  and  Nixon  (1995)  who  investigated  the  role  played  by  turbulence,  and  showed  a  direct  correla¬ 
tion  between  the  shock  motion  and  the  incoming  velocity  fluctuations.  They  also  demonstrated  that 
the  size  of  the  separation  bubble  has,  to  some  extent,  a  weak  effect  on  the  shock  motion.  Despite 
the  lack  of  detail  in  the  inner  layer  (a  log-law  wall  function  was  used  on  a  rough  grid  resolution), 
it  displayed  the  qualitative  features  of  the  shock  oscillation  observed  experimentally  (Dolling  and 
Or  1983). 

A  computation  of  an  adiabatic  turbulent  boundary  layer  flow  past  a  25°  compression  corner  at 
Mach  3.0  and  Res  =  2  x  104  was  performed. 

Allowing  x,  y  and  z  to  denote  the  streamwise,  transverse  and  spanwise  directions,  respectively,  the 
computational  domain  is  Lx  =  16.0(5,  Ly  =  3.4<5,  and  Lz  =  1.925(5.  The  grid  consists  of  213  x  35  x  57 
nodes  in  the  x,  y  and  z  directions,  respectively.  The  reference  quantities  for  non-dimensionalization 
are  length  S  (the  incoming  boundary  layer  thickness),  velocity  U0 0,  density  Poo,  static  temperature 
Too  and  molecular  viscosity  poo  (where  the  subscript  oo  denotes  the  freestream  conditions  upstream 
of  the  compression  corner).  The  tetrahedral  grid  is  employed  and  stretched  in  the  y  direction  with 
a  spacing  of  0.0085  at  the  wall  and  the  stretching  factor  of  1.154.  The  grid  is  concentrated  around 
the  compression  corner.  The  details  of  the  grid  are  shown  in  Table  6,  wherein  A y+  =  Ay/77  with 


Table  6:  Details  of  Grids 


Name 

Mach 

Ax+ 

Ay+ 

at  the  wall 

Az+ 

A  x/S 

A  y/S  A  z/5 

at  y  =  S 

Tetras 

Theoretical  value 

2.88 

24 

1.9 

8.1 

0.1 

0.14 

0.034 

LES 

2.88 

20.9 

1.67 

7.1 

0.1 

0.14 

0.034 

2,018,240 

the  inner  length  scale  77  =  vw/uT  (uw  is  the  kinematic  viscosity  at  the  wall,  uT  =  -*/ tw / pw  is  the 
friction  velocity,  tw  is  the  wall  shear  stress  and  pw  is  the  density  at  the  wall) .  The  theoretical  values 
of  ur  and  vw  are  obtained  from  the  combined  Law  of  the  Wall  and  Wake  evaluated  at  y  =  <5  and 
the  power  law  of  the  relationship  between  temperature  and  kinematic  viscosity,  respectively. 

The  inflow  condition  is  obtained  from  a  separate  flat  plate  boundary  layer  computation.  The  non- 


Figure  23:  Mean  streamwise  velocity  Figure  24:  Mean  temperature 


slip  boundary  condition  is  used  to  the  adiabatic  wall.  All  the  flow  variables  shown  in  the  figures 
are  averaged  in  time  and  the  spanwise  direction.  The  time  averaging  period  is  set  to  three  times 
the  flow-through  time,  where  one  flow-through  time  is  defined  as  the  time  for  the  freestream  flow 
to  traverse  the  computational  domain.  The  details  are  presented  in  Urbin  et  al  1999. 

The  oncoming  flow  characteristics  are  illustrated  by  the  mean  flow  variables  in  Fig.  23  and  Fig.  24 
and  the  Reynolds  shear  stress  in  Fig.  25.  The  comparisons  with  experiments  (Zheltovodov  et  al 
1990)  and  DNS  show  good  agreement. 


Figure  25:  Reynolds  shear  stress 
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Figure  26:  Instantaneous  pressure  contour 

Fig.  26  shows  the  pressure  contour  distribution  at  x  -  y  plane  of  z  =  1.05.  A  strong  separation  and 
attachment  shock  wave  is  formed  at  the  compression  corner  leading  to  the  higher  pressure  level 


after  the  shock.  The  strong  adverse  pressure  gradient  causes  the  skin  friction  coefficient  to  decrease 
dramatically  and  the  flow  separates.  Downstream  of  the  corner,  the  overall  increase  in  pressure  and 
the  decrease  in  Mach  number  cause  the  skin  friction  coefficient  to  recover. 


Figure  27:  Skin  friction  coefficient  Figure  28:  Surface  wall  pressure 


Moo 


Figure  29:  Separation  length  for  LES  and  DNS 

The  computational  results  are  shown  in  Fig.  27-Fig.  29  along  with  experimental  data.  The  skin 
friction  coefficient  in  Fig.  27  is  compared  with  the  experiment  at  higher  Reynolds  number  of  Res  — 
63560.  According  to  the  Law  of  the  Wall  and  Wake,  the  friction  velocity  is  decreased  with  the 
increase  in  Reynolds  number,  leading  to  the  higher  skin  friction  coefficient  in  the  computation.  The 
time  and  spanwise  averaged  surface  pressure  profile  along  the  streamwise  direction  is  compared 
with  experiment  at  higher  Reynolds  number  in  Fig.  28  and  the  pressure  plateau  is  not  observed. 
The  difference  between  the  predicted  and  experimental  surface  pressure  profile  may  be  attributable 
to  the  difference  in  Reynolds  number. 

The  effect  of  Reynolds  number  on  the  separation  length  is  plotted  in  Fig.  29.  In  this  figure,  the 
separation  length  is  measured  by  connecting  the  separation  and  attachment  points  at  which  the 
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time  and  spanwise  averaged  skin  friction  coefficients  go  to  zero  and  then  scaled  by  the  characteristic 
length  ( Lc )  proposed  by  Zheltovodov  and  Schuelein  1987,  1993. 

Lc  =  Sfo/Ppif'/M^  (32) 

where  5  is  the  incoming  boundary  layer  thickness,  p2  is  the  pressure  after  the  shock  in  inviscid  flow, 
ppi  is  the  plateau  pressure  obtained  by  the  empirical  formula  ppi  =  poo(0.5Moo  +  1)  (Zukoski  1967) 
and  Moo  is  the  Mach  number  in  the  uniform  flow.  Some  LES  and  DNS  results  by  other  researchers 
are  also  plotted  in  Fig.  29  for  comparison.  Our  LES  successfully  predicts  the  consistent  trend  with 
the  experiment. 


3  Expansion- Compression  Corner 

Supersonic  expansion-compression  corner  (Fig.  30)  is  reminiscent  of  aerodynamic  configurations 
wherein  a  supersonic  boundary  layer  is  subjected  to  an  initial  expansion  followed  by  a  subsequent 
compression.  Interest  in  this  configuration  is  due  in  part  to  the  stabilizing  influence  of  the  expansion 
(  Dussauge  1987,  Zheltovodov  et  al  1987,  Zheltovodov  and  Schueleinl987,  Smith  1997,  Stephen 
1998,  Zheltovodov  1990).  The  first  systematic  combined  experimental  and  numerical  study  of  an 
expansion-compression  corner  by  Zheltovodov  1992  and  Zheltovodov  1993  showed  that  several 
different  turbulence  models  (including  k—e,  q-ui  and  several  modifications  thereto)  did  not  accurately 
predict  the  separation  and  attachment  positions,  and  distributions  of  surface  skin  friction  and  heat 
transfer.  We  therefore  seek  to  ascertain  the  capability  of  LES  to  predict  this  flowfield. 

An  incoming  Mach  3  adiabatic  equilibrium  turbulent  boundary  layer  of  height  8  expands  over  a 
25°  corner  followed  by  a  25°  compression.  The  distance  along  the  expansion  surface  is  7.15  (i.e., 
the  vertical  distance  between  the  two  horizontal  surfaces  is  35,  and  the  horizontal  distance  between 
the  expansion  and  compression  corners  is  6.435). 

The  Cartesian  coordinates  x,  y  and  z  are  aligned  in  the  incoming  streamwise,  transverse  and  span- 
wise  directions  with  the  origin  at  the  inflow  boundary.  The  computational  domain  is  Lx  =  24.05, 
Ly  =  3.45,  and  Lz  =  1.9255.  The  expansion  corner  is  located  at  45  from  the  inflow  boundary. 
The  grid  consists  of  253  x  35  x  57  nodes  in  the  x,  y  and  z  directions,  respectively,  forming  479,808 
hexahedra  which  axe  subdivided  into  five  tetrahedra  each.  Thus,  the  total  number  of  tetrahedra  is 
2,399,040.  The  grid  is  stretched  in  the  y  direction  with  spacing  0.0085  at  the  wall  and  a  geometric 
stretching  factor  of  1.154.  The  grid  is  concentrated  in  the  streamwise  direction  in  the  neighborhood 
of  the  expansion  and  compression  corners.  The  details  are  shown  in  Table  7  where  A y+  =  AyuT/vw 
where  uw  is  the  computed  kinematic  viscosity  at  the  wall,  uT  =  \frw/pw  is  the  friction  velocity,  tw 
is  the  computed  wall  shear  stress  and  pw  is  the  computed  density  at  the  wall.  The  grid  is  consistent 
with  the  resolution  requirements  for  the  LES  code  established  by  Urbin  2001. 


Table  7:  Details  of  Grid 


Name 

Ax+ 

Ay+ 

at  the  wall 

Az+ 

Ax/8 

Ay/5 
at  y  =  5 

Az/8 

Tetras 

Computed 

20.9 

1.67 

7.1 

0.1 

0.14 

0.034 

2,399,040 

The  inflow  boundary  condition  is  obtained  from  a  separate  flat  plate  boundary  layer  computation. 

Experimental  data  has  been  obtained  by  Zheltovodov  et  al  1987,  Zheltovodov  and  Schuelein  1987, 
Zheltovodov  1990a  and  presented  in  part  in  tabular  form  in  Zheltovodov  1990  for  the  expansion- 


compression  corner  at  Mach  3  and  several  Reynolds  numbers  Res  based  on  the  incoming  boundary 
layer  thickness  5.  The  experimental  conditions  are  listed  in  Table  8,  where  FPBL  and  ECC  imply 
flat  plate  boundary  layer  and  expansion-compression  corner,  respectively.  The  LES  was  performed 
at  a  lower  Reynolds  number  (Res  =  2x  104)  than  the  experiment  (Res  =  4.4  x  104  to  1.94  x  105)  for 
reasons  of  computational  cost.  Additional  LES  cases  will  be  performed  at  higher  Reynolds  numbers. 


Table  8:  Details  of  Experiments  and  Computation 


Cases 

Mach 

Res 

References 

ECC 

2.9 

4.07  x  104 

Zheltovodov  1990 

ECC 

2.9 

6.76  x  104 

Zheltovodov  1990 

ECC 

2.9 

8.0  x  104 

Zheltovodov  1990 

ECC 

2.9 

1.94  x  105 

Zheltovodovl987  and  Zheltovodov  1990a 

ECC 

2.88 

2.0  x  104 

Present  computation 

FPBL 

2.88 

1.33  x  105 

Zheltovodov  1990a 

The  structure  of  the  flowfield  is  shown  in  Figs.  31  and  Fig.  32  which  display  the  mean  static 
pressure  and  streamlines  at  z  —  5.  The  flow  expands  around  the  first  corner,  and  recompresses  at 
the  second  corner  through  a  shock  which  separates  the  boundary  layer  as  evident  in  Fig.  32.  The 
flowfield  structure  is  in  good  agreement  with  the  results  of  Zheltovodov  1987,Zheltovodov  1988, 
Zheltovodov  1990  and  Zheltovodov  1990a  which  are  shown  qualitatively  in  Fig.  30. 
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Figure  31:  Mean  static  pressure  (s  is  separation,  Figure  32:  Mean  streamlines  (s  is  separation,  A 
A  is  attachment)  is  attachment) 

The  mean  velocity  profiles  in  the  x-direction  are  shown  in  Fig.  33  at  x  =  25  and  x  =  65,  where 
x  is  measured  from  the  inflow  along  the  direction  of  the  inflow  freestream  velocity  (Fig.  31).  The 
abscissa  is  the  component  of  velocity  locally  parallel  to  the  wall,  and  the  ordinate  is  the  distance 
measured  normal  to  the  wall.  The  first  profile  is  upstream  of  the  expansion  corner  which  is  located 
at  x  =  45,  and  the  second  is  downstream  of  the  expansion  fan  and  upstream  of  the  separation  point. 
The  computed  mean  velocity  profile  at  the  first  location  is  slightly  fuller  than  the  experiment.  This 
is  consistent  with  the  experimentally  observed  dependence  of  the  exponent  n  in  the  power-law 
U/Uoo  =  (y/5)1/"  on  the  Reynolds  number.  The  second  profile  shows  a  significant  acceleration  of 
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the  flow  in  the  outer  portion  of  the  boundary  layer  due  to  the  expansion. 


Figure  33:  Mean  velocity  Figure  34:  Separation  length 

Zheltovodov  1988  and  Zheltovodov  1993  developed  an  empirical  correlation  for  the  separation  length 
(defined  as  the  minimum  distance  between  the  mean  separation  and  attachment  points  on  the  wall) 
in  the  expansion-compression  corner  interaction.  The  scaled  separation  length  Lsev/Lc  is  observed 
experimentally  to  be  a  function  of  Re$  where  the  characteristic  length  ( Lc )  is  defined  by 

Lc  =  Se(p2/ppl)31/M3e  (33) 

where  6e  is  the  incoming  boundary  layer  thickness  (upstream  of  the  expansion  corner),  p2  is  the 
pressure  after  the  shock  in  inviscid  flow,  ppi  is  the  plateau  pressure  from  the  empirical  formula 
Ppi  —  Pe{\Me  +  1)  where  pe  and  Me  are  the  static  pressure  and  freestream  Mach  number  upstream 
of  the  compression  corner  and  downstream  of  the  expansion  fan.  In  the  computation,  the  location 
is  taken  to  be  x  =  6J.  The  values  of  Me  and  p2  have  been  computed  using  inviscid  theory.  Also, 
Rese  =  1.8  x  104  for  LES  (Rese  =  PeUebe/Pe,  where  pe,[/e  and  pe  are  computed  using  inviscid 
theory).  The  experimental  data  correlation  of  Zheltovodov  1988  and  the  computed  result2  for 
the  scaled  separation  length  is  shown  in  Fig.  34.  The  computed  value  is  consistent  with  a  linear 
extrapolation  of  the  expermental  data. 

The  surface  pressure  profile  in  Fig.  36  displays  a  pressure  plateau  on  the  compression  face  generated 
by  the  separation  bubble.  The  experiments  exhibit  a  trend  of  increase  in  the  size  of  the  pressure 
plateau  region  with  decreasing  Reynolds  number.  The  experimental  data  at  the  lowest  Reynolds 
number  ( Res  =  4.1  x  104)  shows  close  agreement  with  the  computed  results  for  Res  =  2  x  104  for  the 
location,  extent  and  magnitude  of  the  pressure  plateau.  Moreover,  the  shape  of  the  experimental 
pressure  plateau  shows  little  variation  for  Res  <  6.8  x  104,  thus  suggesting  that  the  computed 
pressure  plateau  region  (for  Res  =  2  x  104)  is  accurate.  The  computed  recovery  of  the  surface 
pressure  is  more  rapid  than  in  the  experiment,  however. 

2The  uncertainty  in  the  computed  value  of  Lscp/Lc  is  associated  with  the  uncertainty  in  determining  Se.  We  have 
used  the  streamwise  Reynolds  stress  (<<  p  >><<  v!v!  >>)  to  determine  Se  (Fig.  35),  where  u '  is  the  fluctuating 
velocity  parallel  to  the  wall. 


<p><u,u,>/pcollfi 


Figure  35:  Reynolds  streamwise  stress  Figure  36:  Surface  pressure 


The  computed  and  experimental  mean  skin  fric¬ 
tion  coefficient  c/  =  Tw/^pooU^  are  shown  in 
Fig.  37.  The  computed  separation  and  attach¬ 
ment  points  are  evident.  The  skin  friction  rises 
rapidly  downstream  of  attachment.  The  com¬ 
puted  results  at  Res  =  2  x  104  are  in  close 
agreement  with  the  experimental  data  at  Res  = 
8.0  x  104  and  1.94  x  105  in  the  region  downstream 
of  reattachment. 


Figure  37:  Skin  friction  coefficient 


4  square  Jet 


Turbulent  round  and  plane  jets  are  simple  inhomogeneous  flows  that  can  be  served  to  verify  models 
for  complex  flows  and  have  been  experimentally  and  numerically  studied  extensively  (Panchapake- 
san  1993,  Rodi  1980).  Recently,  noncircular  jets  have  been  gained  much  interest  in  passive  control 
due  to  their  enhanced  jet  mixing  properties  (Grinstein  1992,  Grinstein  1995a,  Grinstein  1995b, 
Grinstein  2001,Gutmark  1999).  A  square  jet  at  a  Reynolds  number  of  3200  and  a  Mach  number 
of  0.3  is  simulated.  Temporal  evolutions  are  visualized  to  characterize  the  dynamics  of  deforming 
vortex  rings,  ribs  and  their  interactions.  Statistical  quantities  are  quantified  and  compared  with 
the  DNS  results  of  Grinstein  et  al  1995. 

The  grid  consists  of  65  x  65  x  65  hexahedral  cells  in  an  unstructured  grid  covering  a  computational 
domain  of  5 D  in  the  streamwise  direction  and  ±3 D  along  the  transverse  directions.  Each  hexahedral 
cell  is  divided  into  five  tetrahedral  cells,  yielding  a  total  of  1.3M  tetrahedra.  A  uniform  grid  is  used 
along  the  streamwise  direction  with  the  hexahedral  grid  spacing  of  Ax/D  =  0.078,  which  is  larger 
than  0.04  used  by  Grinstein  et  al.  in  their  DNS  (Grinstein  1995).  The  grids  are  stretched  along  the 
other  two  directions  and  the  minimum  grid  spacings  are  Ay/D  =  Az/D  =  0.0375.  The  imposed 
boundary  conditions  include  inflow,  outflow  and  wall  boundaries.  At  the  inflow,  the  streamwise 
velocities  are  prescribed  as 

u  =  U[1  +  A  sin(27r  ft)]  (34) 

and 

U  =  0.5  Uoo  [1  -  tanh[M2|y|/D  -  D/(2\y\))]  x  0.5  U, [1  -  tanh[&2(2N/D  -  D/{2\z\))\,  (35) 

where  A  =  0.02  is  the  perturbation  amplitude,  /  is  the  forcing  frequency  (/  =  0.5),  62  =  0.25 Ri/6, 

where  Ri/0  =  40,  R  =  D/2  and  9  is  the  momentum  thickness.  Zero-gradient  condition  is  imposed 
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at  the  outlet  and  symmetry  boundary  conditions  are  used  at  the  side  walls. 

The  iso-surfaces  of  the  total  vorticity  u  =  +  u;|  +  corresponding  to  co  =  0.25upeak  are  shown 

in  Fig.  38.  Azimuthal  nonuniformities  make  the  evolution  of  the  jet  shear  layer  more  complicated 
relative  to  circular  jets.  Close  to  the  jet  exit,  a  smooth  square  vortex  sheet  can  be  observed,  and 
subsequently  rolled-up  vortex-ring  structures  form  due  to  shear-layer  Kelvin-Helmholtz  instability. 
However,  the  vortex  rings  further  downstream  deform  to  non-planar  shape  due  to  self-induction 
mechanism  caused  by  azimuthal  nonuniformities.  The  deformed  vortex  rings  are  connected  with 
the  four  corners  of  the  initial  square  sheet  by  ribs.  The  hairpin  braid  vortices  aligned  with  the 
corners  progress  faster  in  the  diagonal  direction  than  the  others,  “which  results  in  redistribution 
of  energy  between  azimuthal  and  streamwise  vortices”  (Grinstein  2001).  Further  downstream,  the 
jet  development  is  characterized  by  the  strong  interaction  between  vortex  rings  and  braid  vortices, 
which  leads  to  a  final  breakdown  of  the  large-scales  coherent  structures  and  transition  to  the  tur¬ 
bulent  flow.  The  self-induced  deformation  of  the  rings  and  rib  pair  were  explained  to  be  the  leading 
mechanism  for  larger  entrainment  properties  in  non-circular  jets  relative  to  circular  jets  (Grinstein 
1992).  “The  interactions  between  the  streamwise  vortices  and  the  vortex  rings  is  reminiscent  of  the 
interaction  between  ribs  and  spanwise  rollers  in  the  mixing  layer”  (Grinstein  1992).  Mixing  of  jets 
with  surroundings  can  be  enhanced  through  controlling  the  formation,  development  and  interaction 
of  large-scale  coherent  structures  passively. 

Fig.  39  shows  the  instantaneous  crosswise  vorticity  uz  =  dv/dx  -  du/dy  contours  at  the  central 
x  —  y  plane.  Rolled-up  structures  can  be  observed  near  the  base,  and  subsequently  symmetrical 
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counterrotating  toroidal  structures  form  in  the  shear  layer  which  are  then  followed  by  their  stretch¬ 
ing  and  deformation.  The  organized  structures  can  be  broken  down  into  smaller  eddies  further 
downstream.  Evidently,  the  MILES  model  can  capture  the  transition  process.  Large-scale  vortex 
rings  dominate  in  the  near-field  and  the  small  vortices  dominate  downstream  after  the  breakdown. 
The  spatial  spreading  can  be  clearly  observed  downstream  as  the  jet  spreads  by  entraining  mass 
from  the  surrounding  nonvortical  fluid. 

Fig.  40  (a-d)  shows  the  contours  of  instantaneous  streamwise  vorticity,  ljx  =  dv/dz  -  dw/dy,  across 
the  y  —  z  planes  of  x/D  =  1,  2,  3  and  4.  Quite  different  behaviour  can  be  observed  at  different  axial 
positions.  Vortex  shears  with  some  rounded-corners  are  stretched  and  thickened  but  still  keep  the 
initial  square  shape  at  the  position  of  x/D  =  1.  The  jet  cross  section  switches  axis  45°  relative  to 
that  of  the  jet  nozzle  at  the  axial  location  of  x/D  =  2, 3  due  to  self-induced  velocity  around  the 
corners  by  the  presence  of  streamwise  vorticity.  The  flow  structure  develops  into  a  irregular  shape 
further  downstream  at  x/D  =  4. 

The  statistical  quantities  are  obtained  by  averaging  over  five  forcing  cycles  after  a  statistically 
stationary  state  has  been  reached  after  an  elapsed  dimensionless  physical  time  of  ten.  The  centerline 
distributions  of  the  mean  axial  velocity  by  LES  compared  with  the  DNS  results  of  Grinstein  et  al 
1995  are  plotted  in  Fig.  41.  The  mean  velocity  initially  decays  within  the  first  1.5  diameters  and 
subsequently  shows  a  slight  increase  due  to  the  periodic  roll-up  by  the  sinusoidal  forcing.  The 
decay  after  3.2  diameters  is  the  result  of  turbulent  mixing.  In  general,  the  agreement  between 
present  MILES  and  previous  DNS  results  are  good. 

The  corresponding  centerline  r.m.s.  velocities  v! jUc  are  shown  in  Fig.  42.  Note  that  close  to  the 
inflow  plane  the  r.m.s.  of  velocity  fluctuations  is  about  2  percent  and  corresponds  to  the  imposed 
disturbance  level.  In  the  potential  core  region,  the  r.m.s.  velocity  decreases  slightly  with  increasing 
axial  distance  which  shows  a  tendency  to  remain  laminar  with  low  turbulence  intensities.  Beyond 
the  end  of  the  potential  core,  the  fluctuations  of  velocity  increase  very  rapidly,  which  indicates 
the  appearance  of  the  secondary-instability  mechanism  which  leads  to  the  final  breakdown  of  the 
large  vortex  structures.  Due  to  insufficient  length  scale  in  the  streamwise  direction,  the  self-similar 
behaviour  has  not  been  achieved.  However,  it  can  be  seen  that  the  agreement  between  MILES  and 
DNS  is  very  good  for  the  transient  square  jet. 
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